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ABSTRACT 

We present time series photometry of the M dwarf transiting exoplanet system 
GJ 436 obtained with the EPOCh (Extrasolar Planet Observation and Characteri- 
zation) component of the NASA EPOXI mission. We conduct a search of the high- 
precision time series for additional planets around GJ 436, which could be revealed 
either directly through their photometric transits, or indirectly through the variations 
these second planets induce on the transits of the previously known planet. In the case 
of GJ 436, the presence of a second planet is perhaps indicated by the residual orbital 
eccentricity of the known hot Neptune companion. We find no candidate transits with 
significance higher than our detection limit. From Monte Carlo tests of the time se- 
ries, we rule out transiting planets larger than 1.5 interior to GJ 436b with 95% 
confidence, and larger than 1.25 with 80% confidence. Assuming coplanarity of ad- 
ditional planets with the orbit of GJ 436b, we cannot expect that putative planets with 
orbital periods longer than about 3.4 days will transit. However, if such a planet were 
to transit, we rule out planets larger than 2.0 R^ with orbital periods less than 8.5 days 
with 95% confidence. We also place dynamical constraints on additional bodies in the 
GJ 436 system, independent of radial velocity measurements. Our analysis should serve 
as a useful guide for similar analyses of transiting exoplanets for which radial velocity 
measurements are not available, such as those discovered by the Kepler mission. From 
the lack of observed secular perturbations, we set upper limits on the mass of a second 
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planet as small as 10 in coplanar orbits and 1 in non-coplanar orbits close to 
GJ 436b. We present refined estimates of the system parameters for GJ 436. We find 
P = 2.64389579 ± 0.00000080 d, = 0.437 ± 0.016 Rq, and Rp = 3.880 ± 0.147 iZ®. 
We also report a sinusoidal modulation in the GJ 436 light curve that we attribute to 
star spots. This signal is best fit by a period of 9.01 days, although the duration of the 
EPOCh observations may not have been long enough to resolve the full rotation period 
of the star. 

Subject headings: eclipses — stars: individual (GJ 436) — stars: planetary systems — 
techniques: image processing — techniques: photometric 



Introduction 



EPOXI (EPOCh + DIXI) is a NASA D iscovery Program Mission of Opportunity using the 
Deep Impact flyby spacecraft (iBlumd l2005l ). From January through August 2008, the EPOCh 
(Extrasolar Planet Ob servation and Characterization) Science Investigation used the HRI camera 
(jHampton et al.ll2005l ) with a broad visible bandpass to gather precise, rapid cadence photometric 
time series of known transiting exoplanet systems. The majority of these targets were each observed 
nearly continuously for several weeks at a time. In Table 1 we give basic information about the 
seven EPOCh targets and the number of transits of each that EPOCh observed. 

One of the EPOCh science goals is a search for additional planets in these systems. Such 
planets would be revealed either through the variations they induce on the transits of the known 
exoplanet, or directly through the transit of the second planet itself. This search is especially 
interesting in the case of the GJ 436 systern , since the non-zero eccentricity of the known Neptune- 
mass planet, first measured bvlButler et al.l (|2004l ). may indicate the presence of a second planetary 
companion (jManess et al.l 120071 ). Because GJ 436 is a nearby M dwarf, it is also the only EPOCh 
target for which we are sensitive to planets as small as 1.0 We will describe the searches for 
additional planets conducted on the remaining EPOCh targets in subsequent papers. 

The search for transiting Earth-sized planets in the GJ 436 light curve is scientifically com- 
pelling for the following four reasons. First, theoretical predictions of the mass-radius relation for 
"super Earths" are still largely observationall y unconstrained, with the excitii ig exceptions of the 
two k nown transiting super Earths CoRoT-7b (ILeger et al.lboogi ') and GJ 1214b (ICharbonneau et al 



20091 ). Depending on the level of observational uncertainty, a measurement of the mass and radius 
of a super Earth could point to the presence of a large amount of water or iron (enabled with 10% 
uncertainty), or allow us to distinguish betwee n a planet compos ed predominately of water ice, 
silicates, or iron (enabled with 5% uncertainty; ISeaeer et al.l 120071 ). Second, the discovery of two 
transiting bodies in the same system would permit the direct observation of their mutual dynami- 
cal interactions. This would enable constaints on the masses of the two bodies independent of any 
radial velocity measurement (jHolman &: Murravll2005l : lAgol et al.ll2005l ). Since radial velocities can 
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only be observed for planets above a certain mass limit, this is an important tool for future surveys 
of stars too faint for radial velocity measurements. Third, the discovery of an Earth-sized planet 
at an orbital radius outside that of a giant planet would inform theories of planet formation. Hot 
Earths are predicted to be captured in low order mean motion resonances with migrating giant 



planets (IRavmond et al 



2006 



Mandell et al.l 120071 ). Since the phenomenon of Earth-sized planets 



at larger orbital radii than Jovian planets is not observed in our own solar system, observations of 
exoplanet systems are particularly important for this question. 



Finally, the eccentricity of the known transiting Neptune-mass planet, GJ 436b (jButler et al 



2004 ) ■ may indicate the presence of an additional perturbing planet, since the assumed circulariza- 
tion timescale for the known planet is m u ch less than the ag e of the system (jManess et al.l 120071 : 



Deming et al.l 120071 : iDemorv et al.l 120071 ) . iRibas et al.l (|2008l ) claimed evidence for a 5 super 



Earth in radial velocit y observations of GJ 436, but this p roposed planet was ruled out by sub 



sequent investigations ( Alonso et al. 2008; Bean Seifahrtl[2008 ). The absence of this additional 



perturbing body in the GJ 436 system would also be very scientifically interesting. If no other 
body is present to explain the eccentricity of GJ 436b, the observed eccentricity requires a very 
high tidal dissipation parameter, Q. The current estimate of the circularization timescale assumes 
a Q value for the hot Neptune similar to the value derived for the ice giant planets in our own solar 
system, so a substanti ally different Q wou ld indicate physical pro perties of GJ 436b very different 
from these ice giants (jPeming et al.l 120071 ). iJackson et al.l (120081 ) show that a ratio of planetary 
tidal dissipation parameter to planetary Love number Q/k2 for GJ 436b greater than 10^'^ can 
explain the system's eccentricity (the Love number /c2 is theoretically between 3/2, in the case of a 
homogeneous body, and 0, in the c ase of a cent rally condensed body, but ranges between 0.3 and 0.6 
for gas giants in the Solar Svstem: lBursalll992l l. In contrast, Uranus and Neptune, the solar system 
bodies presumably most similar in composition and mass to GJ 436b, have tidal Q parameters esti- 
mated at Qtj < 3.9 X 10"^ and 1.2 x 10^ < Qtv < 3.3 x 10^ respectively dTittemore fc Wisdoml[l989l : 
Banfield &: Murravlll992l ) — several orders of magnitude smaller than the Q necessary to explain 
the eccentricity of GJ 436b . If the eccentricity is not attributable to a high Q, there may instead be 
an additional perturbing body maintaining the system's eccentricity. The possibility of a close-in 
resonant comp anion in 2:1 or 3:1 r esonance with G J 436 b is strongly disfavored by transit timing 
measurements (jPont et al.l l2009l ) . iBatvgin et al.l (j2009l ) proposed possible secular perturbers to 



GJ 436b, the presence of which would be consistent with o bserved radial velocities , transit timing 
measurements, and the non-zero eccentricity of the system. iBean fc SeifahrtI (j2008l ) also quantified 
the improvement to the goodness-of-fit of the GJ 436 radial velocity data with the addition of 
perturbing planets to the model — the parameter space t hey investigated inclu ded putative planets 
of lower mass and eccentricity than those suggested by iBatvgin et al.l (j2009l ) . The existence and 
possible orbital parameters of this putative planet are still open questions. In addition, the recent 
discovery of the second transiting ho t Neptune, HAT-P -llb, also makes this question timely, since 
the planetary orbit is also eccentric (jBakos et al.ll2010l ). 



The remainder of this paper is organized as follows. In Section 2, we describe the photometry 
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pipeline we created to produce the time series. In Section 3, we detail the refinement of system 
parameters and the search we conduct for additional planets around GJ 436, both for additional 
transits and for dynamical perturbations to GJ 436b. We present a Monte Carlo analysis of the 
EPOCh observations of GJ 436 and demonstrate the sensitivity to detect a transiting planet as 
small as 2.0 times the size of Earth for all periods less than 8.5 days with high confidence. We 
discuss the upper limits on the mass of additional coplanar and non-coplanar planets with periods 
between 0.5 and 9 days from dynamical constraints. We also discuss the constraints we place on the 
rotation period of GJ 436. In Section 4, we present our best candidate transit signal, and from the 
search for additional transits we place upper limits on the radius and mass of the putative planet 
GJ 436c. 



2. Observations and Data Reduction 



We acquired observations of GJ 436 nearly continuously during 2008 May 5-29, interrupted 
for several hours at approximately 2-day intervals for data downloads. The basic characteristics 
of the targets and observations are given in Tables 1 and 2. Observations of this type were not 
contemplated during development of the original Deep Impact mission; the spacecraft was not de- 
signed to maintain very precise pointing over the timescale of a transit (Table [2]). Furthermore, the 
available onboard memory precludes storing the requisite number of full- frame images (1024x1024 
pixels). Hence the observing strategy used 256x256 sub-array mode for those times spanning the 
transit, and 128x128 otherwise. This strategy assured complete coverage at transit, with minimal 
losses due to pointing jitter exceeding the 128x128 sub-array at other times. 

We use the existing Deep Impac t data reduction pipe line to perform bias and dark subtractions. 



as well as preliminary flat fielding (jKlaasen et al.ll2005l ). The image motion from pointing jitter 



produces a significant challenge for photometry at our desired level of precision. The flat field 
calibration that was obtained on the ground before launch is not successful at the level of precision 
needed here, because spatial variation of the sensitivity of the CCD has changed modestly since 



Table 1. EPOCh Targets 



Name 


V Magnitude 


Number of Transits Observed^ 


Dates Observed [2008] 


HAT-P-4 


11.22 


10 


Jan 22-Feb 12, Jun 29-Jul 7 


TrES-3 


11.18 


7 


Mar 8-March 10, March 12-Mar 18 


XO-2 


12.40 


3 


Mar 11, Mar 23-Mar 28 


GJ 436 


10.67 


8 


May 5-May 29 


TrES-2 


11.41 


9 


Jun 27-Jun 28, Jul 9-Jul 18, Jul 21-Aug 1 


WASP-3 


10.64 


8 


Jul 18-Jul 19, Aug 1-Aug 9, Aug 11-Aug 17 


HAT-P-7 


10.50 


8 


Aug 9-Aug 10, Aug 18-Aug 31 



^Some transits are partial. 
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launch in 2005. Our observing sequences included observations of a green stimulator LED ("stim 
lamp") that could be switched on to illuminate the CCD. The stim lamp illumination is neither 
flat nor stable in an absolute sense, but its spatial pattern was designed to be stable. Hence it is 
useful to define and monitor changes in the flat field response pattern of the CCD. But since the 
stim lamp has a different effective wavelength than the stars, it is not a perfect calibrator for flat 
field changes. For this reason we also use a 2D spatial-spline fit to the actual data, as a bootstrap 
fiat field method as described below. 

We extract the photometric time series as follows. We determine the position of the star 
on the CCD using PSF fitting, by maximizing the goodness-of-fit (with the statistic as an 
estimator) between an image and a model PSF with variable position, additive sky background, 
and multiplicative brightness scale factor. We take advantage of a defocused aperture; given our 
limited on-board memory, defocusing enables us to spread the starlight over more pixels and extend 
our duty cycle. The PSF itself has a donut-like shape with a roughly 10 pixel FWHM. A model of 
the PSF is produced from the drizzle of more than 1200 60x60 pixel cutouts, filtered to eliminate 
cosmic ray hits before drizzle. The final PSF model is sampled to a tenth of a pixel. A bilinear 
interpolation of this PSF increases the sampling to a hundredth of a pixel, which is the accuracy 
to which we estimate the position from the grid. At this point, we perform cosmic ray filtering 
by removing images from the sample with a larger than 15a outlier in the residuals between image 
and best-fit PSF model. Because of the high cadence of the EPOCh observations, we simply reject 
the approximately 30 images (RiO.1% of the total) containing a cosmic ray overlying the stellar PSF 
from the time series. 

From each image we subtract a bias determined from the sigma-clipped median of the over- 
clock pixels in each of the four quadrants of the CCD. These are not true pixels that lie on an 

unilluminated part of the CCD, but rather the bias values read out after the true pixels and then 
recorded to pixels on the outside of the FITS images. We subtract a bias independently for each 
quadrant, since the original Deep Impact reduction pipeline does not account for time-dependent 
bias variations of the CCD. 

We then process the images to remove several sources of systematic error. 

Table 2. Characteristics of the EPOCh Observations 



Instrument Parameter 



Value 



Telescope aperture 
Spacecraft memory 
Bandpass 
Integration time 
Pointing jitter 
Defocus 
Pixel scale 
Subarray size 



30 cm 
300 Mb 
350-1000 nm 
50 seconds 
± 20 axc-sec per hour 

4 arc-sec FWHM 
0.4 arc-sec per pixel 
256x256 pixels spanning transit, 128x128 otherwise 
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1. We scale down the two central rows by a constant value. Due to the CCD read-out 

electronics, there is a reduction in signal in the pixels at the internal boundary of the two upper 
and lower imaging regions. However, because the Deep Impact pipeline flat flelds these pixels in 
the same way as the others, we observe the brightness of the star to increase by about 3% when 
the stellar PSF lies on the central rows if we do not apply this correction. 

2. We scale down the central columns by a separate constant value. We observe an increase in 
brightness on the order of 0.25% when the stellar PSF overlies the central columns if this correction 
is not applied. We interpret the physical origin of this sensitivity variation to be the serial read-out 
register, which is split in the middle of the CCD to allow the rows to be read out at both ends of 
the register. 

3. We scale the entire image by a multiplicative factor determined by the size of the sub-array. 
We gathered observations using two sub-arrays of the 1024x1024 CCD as mentioned above: one 
in 128x128 mode and another in 256x256. Wc observed an offset in the average out-of-transit 
brightness between the two sub-arrays of 8.5 x lO'*. We correct for the offset in 256x256 mode by 
performing the photometric extraction of the time series, determining the decrement in brightness 
observed in 256x256 mode, uniformly dividing the 256x256 images by this value, and repeating 
this process until the out-of-transit brightness shows no offset between the two observing modes. 

4. We divide the images by a "stim" , described above. We first bias-correct the stims using the 
same prescription as we apply to the images. We then process the stims to remove the asymmetrical 
illumination pattern, which we model as a plane surface in x and y position. 

We then perform aperture photometry on the corrected images. We choose an optimal aperture 
radius based on analysis of the standard deviation of the out-of-transit time series. We find that 
this standard deviation is minimized for an aperture radius of 10 pixels, corresponding to twice the 
HWHM of the PSF. After performing the bias subtraction, scaling of the two middle columns and 
rows, scaling of the 256x256 sub-array images, and division by a stim, the time series still suffers 
from significant red noise due to the interpixel sensitivity variations on the CCD. At this point, we 
implement a 2D spline fit to the data with position by fitting a surface, with the same resolution as 
the CCD, to the brightness variations on the array. We randomly draw a subset of several thousand 
out-of-transit and out-of-eclipse points from the light curve (from a data set of ^^29,988 points) 
and find a robust mean of the brightness of the 30 nearest neighbors for each. We then fit a spline 
surface to these samples, and correct each data point individually by linearly interpolating on this 
best-fit surface. We use only a small fraction of the observations to create the spline surface in 
order to minimize the potential transit signal suppression introduced by flat flelding the data by 
itself. 

To produce the final time series, we iterate the above steps, fitting for the row and column 

multiplicative factors, the 256 mode scaling factor, and the 2D spline surface that minimize the 
out-of-transit white noise of the photometric time series. We include one additional step to create 
the final 2D spline, which is to iteratively remove an overall modulation from the GJ 436 light 
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curve, which has a roughly sinusoidal shape with an amplitude of a few parts in 10^. We attribute 
this modulation to star spots, and discuss this signal in Section 3. After performing the 2D spline, 
we fit a polynomial to the corrected and binned flux, divide this polynomial from the pre-splined 
time series, and repeat the spline fit. Otherwise, we expect the modulation signal to introduce 
red noise to the time series, since we correct for interpixel variations with the assumption that the 
star's intrinsic brightness is constant outside of times of eclipse and transit. 

After we take these steps to address the systematics associated with the observations, the red 
noise is largely removed. Figure [1] shows the GJ 436 time series before and after this 2D spline 
correction. In the bottom panel, we show that the time series after the 2D spline bins down roughly 
as predicted for Gaussian noise over timescales less than 4 hours. In the corrected time series, we 
attribute the scatter to photon noise and low- level cosmic rays. We compare the sigma-clipped 
standard deviation of the out-of-transit and out-of-eclipse flux to the expected value of the photon 
noise-limited precision, and find that we are 56% above the Poisson limit. 

~-\ 1 r 1 1— — n 1 1—: — i — — i r — r 1 — n 1 r": — i 1 r 1 1 — 

: EPOCh GJ 436: ■ . ■ . ■ I ■ ■ ' - , . After : _ 

: I : i : [ : j : I : 1 : i : i : ; 

j_j !_i L L! I J ii I : I Li L I I J I -Lj I I I t I I !_!_! 

5 10 15 20 

Time [Days after First Transit Prior to Observations] 

Bin Size [Hours] 

0.1 1.0 10.0 100.0 

— I — I — I — I I I 1 1 1 1 — I — I — I I I 1 1 1 1 — I — I — I I I 1 1 1 1 — I — I — I I I 1 1 1 — : 
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Fig. 1. — Top panel: GJ 436 time series before (lower curve) and after (upper curve) 2D spline 
correction. The uncorrected time series (lower curve) has had the two middle rows and columns in 
each image, and the entire image if taken in 256x256 observing mode, scaled by a multiplicative 
factor to reduce the fiux dependence on position and sub-array size. The images were also divided 
by a flat fleld constructed from a stim. We have used the 2D spline to correct for additional 
interpixel variation in the upper curve. Bottom panel: The data (diamond symbols) bin down 
consistently with the expectation for Gaussian noise (shown with a line, normalized to match the 
value at N=l). 

We also investigate the transit signal suppression introduced by using a flat field created from 
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the out-of-transit and out-of-eclipse data itself. We avoid the suppression of transits of GJ 436b 
by excluding those observations from the points used to generate the flat field surface, so that we 
only use the presumably constant out-of-transit and out-of-eclipse observations to sample the CCD 
sensitivity. However, if the transit of an additional planet occurs while the stellar PSF is lying 
on a part of the CCD that is never visited afterward, the 2D spline algorithm instead treats the 
transit as a dark pixel. To quantify the suppression of additional transits that result from using the 
2D spline, we inject transit light curves with periods ranging from 0.5 days to 20 days in intervals 
of 30 minutes in phase (ranging from a phase of zero to a phase equal to the period) into the 
EPOCh light curve just prior to the 2D spline step. After performing the 2D spline, we then phase 
the data at the known injected period and fit for the best radius, using as the goodness-of-fit 
statistic. We find that the best-fit radius is suppressed to a mean value of 73% at all periods, with 
the standard deviation from that value increasing with period from 3% at at period of 1 day to 
16% at a period of 10 days. We describe our incorporation of signal suppression into our search for 
additional planets in greater detail in Section 3.2. 



3. Analysis 

3.1. Transit Parameters 

We describe here our refinement of the GJ 436 system parameters. When conducting our 
analysis, we were careful to account for the effects of remaining correlated noise on the parameters 
and their uncertainties. In the final calibrated GJ 436 light curve, we still observe evidence of 
correlated noise and trends that have not been corrected by the reduction process. For each 
transit, we fit a line to the out-of-transit data on both sides of the transit (from 3 minutes outside 
of transit to half an hour outside of transit) and divide the time series by this line. 

We investigate the effects of limb-darkening using several different methods. In the first in- 
stance, we use stellar atmosphere models to fix the limb-darkening coeffi cients to a set of the oretical 
values. Initially, we use a model atmosphere produced by R. Kurucz (lKurucall994l . |2005| ) corre- 
sponding to Teff = 3500K, log g = 4.5, [M/H] = 0.0 and Vtur h = 0.0 km/s. We fit the four 
coefficients of the non-linear limb-darkening law of IClaretl (j200d ) to 17 positions across the stellar 
limb. We repeat this fit in 0.2 nm intervals across the EPOXI bandpass, weighted for the total 
sensitivity (including filter, optics and CCD response) and photon count at each wavelength. We 
calculate the final set of coefficients as the average of the weighted sum across the bandpass, for 
which we find ci = 0.97, C2 = —0.50, C3 = 0.54, and C4 = —0.17. In order to understand the 
effect of the stellar atmosphere model choice on the fi nal derived transit para meters, we also gen- 
erate a set of coefficients from the PHOENIX model (|Hauschildt et al.lll99d ) with Teg = 3300i^, 
log g = 4.5, [M/H] = 0.0 and i^turb = 0.0 km/s (private communication), and find ci = 0.97, 
C2 = 0.35, C3 = —0.21, and C4 = —0.03. As a final test, we also run a parallel analysis where 
quadratic limb-darkening coefficients are allowed to vary in the fit described below, to see if the 
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data are sufficiently precise to break the degeneracy of the transit shape with geometric parameters 
and hmb-darkening. This is discussed further below. 

For each of the three limb-darkening treatments, we fit the transit using the analytic algorithms 
of lMandel &: Agoll ( 20021 ) . Using the Levenberg-Marquardt algorithm with as the goodness-of- 
fit estimator, we fit for t h ree ge ometric parameters: Rp/Ri,, Ri,/a and cos i. Initially we fix the 
period to the iBean et al.l (|2008l ) value of P = 2.643904 days and allow the times of the center of 
transit to vary independently for each of the eight transits, for a total of eleven free parameters. 
In addition, there are two quantities that parametrize quadratic limb darkening. These are set as 
two orthogonal linear combinations of the coefficients 271 + 72 and 272 — 71 • In all three cases we 
find results that are internally consistent at the la level, demonstrating that the choice of limb- 
darkening treatment does not significantly alter the derived parameters. We note that we can only 
constrain one linear combination of the limb-darkening coefficients in the quadratic limb-darkening 
case, the other being degenerate with the geometric parameters to within the precision of the light 
curve. Due to this degeneracy, the error bars on the geometric parameters are larger in the case 
where the limb-darkening coefficients are allowed to vary; the uncertainty on the planetary radius 
Rp is larger by a factor of 2.5 and the uncertainty on the stellar radius i?^ is larger by a factor of 
2.2. 

For the Kurucz stellar atmosphere limb-darkening coefficients, we find Rp/R^ = 0.08142 =b 
0.00085, R^/a = 0.0707 ±0.0025 and cos i = 0.0594 ±0.0030. The light curve is phased and binned 
into 2 minute intervals with this best-fit model overplotted and shown in Figure [2l Since there 
are significant time-correlated systematic errors in the light curve, we c alculate the errors on the 
parameters using the "rosary bead" method, in the fashion described bv lWinn et al.l (j2008l ). First, 
we subtract the best-fit model from the light curve and compile a set of initial residuals. We then 
shift the residuals along the light curve to the next time stamp in each case, preserving the temporal 
correlation, add the model back to the residuals, and fit the new light curve as above. For GJ 436, 
we repeat this process 200 times, which corresponded to a total shift of over three hours, more 
than three times the duration of the GJ 436 transit. This is sufficient to sample the systematics, 
which on short timescales vary on the order of 10-30 minutes. The error bars on the parameters 
are then calculated as the range required to encompass 68% of the results. Increasing the number 
of shifts to the residuals does not increase the resulting error bars. We also use the Monte Carlo 
Markov Chain method, as adapted to transit light curve analysis by e.g. iHolman et al.l (|2006l ). and 
find that it results in significantly smaller error bars (30-40%) than the rosary bead analysis, due 
to the inability to factor in systematic errors in the light curve. The Monte Carlo Markov Chain 
method remains a useful tool, however, for assessing the impact of systematic errors on the derived 
parameters. 



Using the stellar mass from lTorresI ((20071), = 0.452 ±0.013, we find R^ = 0.437 ±0.016 R^ 



'0; 



Rp = 3.880±0.147 i?0, and i = 86.60±0.1 7° . Our mea s uremen t of the plan e tary r adius is consistent 



at the Icr level with the measurements of lCillon et al.l (j2007bl ) . iPont et al.l (|2009l ). and ISouthworth 
(j2008l ). and only marginally inconsistent (at the level of 1.2>a and lAa) with the values measured by 
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Shporer et al. ( 2009) and Gillon et alj ( 2007a). It is sign i ficant ly smaller than the values obtained 
by iTorresI (120071 ) . lOeming et al.1 ((20071), and ISean et al.1 ((20081). Although different treatments of 
limb darkening may partly account for the discrepancy, in all three cases of limb darkening analysis 
described above, we find internally consistent results. Our sm aller value of th e plan etary radius 
would require a reduced mass of the H/He envelope of GJ 436b. ICouehlin et al.l (j2008l ) tentatively 
suggested that the GJ 436 inclination, transit width and transit depth were increasing with time as 
the eccentric orbit precessed. Although we do find an slight increase in inclination, again consistent 
with previous published values, we find a shorter (58.5±2.5 min) and shallower (6.173ib0.037 mmag) 
transit than expected from the predicted increase. These findings weaken the trends observed by 



Coughlin et al.l (j2008l ) and we therefore cannot confirm any parameter variation with these data. 



Phased and Binned GJ 436 Light Curve 




0.01 0.00 0.01 J 



0.01 0.00 0.01 

Phose 



0) 



0.002 
0.001 
0.000 

-0.001 
-0.002 



Fig. 2. — Phased and binned GJ 436 light curve with the best-fit model overplotted. The error 
bars are the rms errors in the bins. 



We fix the times of the centers of transit at the values returned by the least-squares fit de- 
scribed above. Th ese times are shown in Table |3| i n UTC . Adding our e i ght tr a nsit times to 
those pub hshed bv IPont et al.l (l2009l). ICaceres et all hood). Ishporer et al.1 kood). I Gillon et al 



(|2007af lbl). ICoughlin et al.l (|2008l ) . iBean et al.l (|2008l ) . iDeming et all (|2007l ) and lAlonso et al.l (|2008l ) 
we find a new hnear, weighted ephemeris of Tc(BJD) = 2, 454, 600.693205 ± 0.000040 and P = 
2.64389579 ± 0.00000080 d. 
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3.2. Search for Additional Transiting Planets 

We perform a robust search of the GJ 436 time series for evidence of additional planets using 
two methods. The first is a search for additional shallow transits in the light curve. We developed 
software to search f or these additiona l transits partly modeled on the methods employed on the 
MOST photometry (jCroll et al.lbnOTal lbl'l. The steps involved in the procedure are described in this 
section. 

The photometric precision of the GJ 436 time series (with S/N~65 for each transit event) 
should enable a detection of a 1.5 R® planet with S/NwlO and a 2 with S/NsslT, even if the 
planet produces only a single transit event. In order to test this prediction, we conduct a Monte 
Carlo analysis to assess how accurately we could recover an injected planetary signal in the GJ 436 
light curve. We evaluate our sensitivity to transit signals on a grid in radius and period space 
sampled at regular intervals in Rp and regular frequency spacing in P. We create an optimally 
spaced grid as follows: for the lowest period at each radius, we determine the radii at which to 
evaluate the adjacent periods by solving for the radius at which we achieve equivalent signal-to- 
noise (for t his reason, we expect s ignificance contours to roughly coincide with the grid spacing). 



We use the iMandel Sz Agoll (j2002l ) routines for generating limb-darkened light curves given these 
parameters to compute a grid of models corresponding to additional possible planets. If we make 
the simplifying assumptions of negligible limb darkening of the host star, a circular orbit, and an 
orbital inclination angle i of 90°, the set of light curves for additional transiting bodies is a three 
parameter family. These parameters are radius of the planet Rp, orbital period of the planet P, 
and orbital phase (j). At each test radius and period, we inject planetary signals with 75 randomly 
assigned phases into the residuals of GJ 436 EPOCh light curve with the best GJ 436b transit 
model divided out, and then attempt to recover blindly the injected signal by minimizing the 
statistic. We first conduct a coarse gi'id search in radius, period, and phase. We select the 
spacing of this grid to minimize processing time while ensuring that the transit was not missed. 
We sample the space at 300 points in period space (at even frequency intervals between 0.5 and 
8.5 days), 50 points in radius space (between 0.5 and 5.5 Earth radii) and a variable number of 
points in phase space from around 30 to nearly 200 (sampled at half a transit duration for eac h 



period). We use an expression for the transit duration r given bv lSeager &: Mallen-Qrnelaa (|2003l ): 



sm I sm 



/vrr 




cos^i. 



(1) 



For each test model, we compute the using the out-of-transit standard deviation to estimate 
the error in each point. The grid search requires about 24 CPU hours on a 2.66 GHz processor 
for each radius and period (75 light curves with randomly inject ed phases). After th e grid 



minimum is determined, we use the amoeba minimization routine (jNelder &: Meadlll965l ) to more 
finely sample the x"^ space in order to find the x^ minimum from the specified nearest grid point. 
We also investigate whether aliases of the best-fit period from the x^ g^id improve the fit. We 



- 12 - 



find that roughly half of the best solutions from the grid are aliases of the injected period, most at 
either half or twice the value of the injected period, but we test aliases at every integer ratio from 
1/35 to 35 times the given period (these other aliases occur less than one percent of the time). We 
also repeat the finer grid search at the three next lowest minima, in case the best solution (or 
an alias of the best solution) lies closer to that grid point. For all 3600 injected signals, we recover 
a solution which is a better fit (in the sense) than the exact injected signal. For injected radii 
greater than 1 R^, we recover the period to within 1% for all injected periods less than 7 days 
(with the exception of one signaQ). For these reasons, we are confident that we are sampling the 
space sufficiently finely to locate the best solution. 

We quantify the success of this analysis by how well the search blindly recovers the known 
injected transit signal. We define the error on the recovered parameter, for instance period, to be 
I Pinjected — Pobserved \ / Pinjected- We Set an approximate error of one sigma for a given parameter 
at the value that includes 51 of the 75 errors at that point in radius and period (roughly 68% of 
all values). Figure [3] shows this relative error in radius for all searches. We shift the location of 
the points at which we evaluate our sensitivity, shown as diamonds, to conservatively incorporate 
the expected level of signal suppression from a planetary transit. As we note in the last paragraph 
of Section 1, we anticipate suppression of additional transit signals from the bootstrap fiat field 
treatment of the EPOXI data. We evaluate the suppression we expect at all eight periods on the 
grid shown in Figure El We inject planetary transits into the pre-flat fielded light curve in intervals 
of 30 minutes in phase (from a phase of zero to a phase equal to the period) , apply the 2D spline 
flat fleld, phase to the known period, and fit for the suppressed transit depth. In this way, we 
obtain a complete understanding of the suppression at each period, and we can evaluate the level 
of suppression that we expect with a given confidence. In particular, with 95% confidence, we 
expect suppression to be no worse than 71% the injected radius at 0.55 days, 68% at 0.79 days, 
65% at 1.13 days, 67% at 1.63 days, 63% at 2.35 days, 65% at 3.38 days, 53% at 4.86 days, and 
55% at 6.98 days. We incorporate this expected suppression by shifting the effective radius values 
of the grid points at which we evaluate our sensitivity to additional transits. For example, at 1.63 
days, all grid points have been shifted upward in radius by a value of 1/0.67, or 1.49. Because we 
anticipate no more than 67% radius suppression at this period (with 95% confidence), an authentic 
transit signal would pessimistically appear only 67% its original radius value after we process the 
data. For this reason, the recovery statistics corresponding to a 1.0 i?® transit depth in the final 
light curve would be accurate for an original transit signal of a 1.49 planet once we fold our 
expectation of signal suppression. For a planet with 2.0 times the Earth's radius, we are able to 
recover the radius to better than 5% accuracy at least 68% of the time for all periods less than 4 
days, and better than 10% at least 95% of the time in the same period interval. 



In this injected signal with a period of 7 days, only one transit event occurred over the baseline of the observations. 
Although the radius of the planet corresponding to the single transit event was recovered to 1%, more than a single 
transit event is required to accurately recover the period itself. 
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Fig. 3. — Constraints on radius from the Monte Carlo analysis. For each point in radius and period, 
we create 75 light curves with random orbital phases, inject them into the GJ 436 residuals, and 
attempt to recover them blindly. The diamonds indicate the grid of radii and periods at which we 
evaluate our sensitivity; the contours are produced by interpolating between these points. At left, 
contours indicate the relative error in radius (absolute value of recovered-injected/injected radius) 
that encloses 68% of the results. At right, contours indicate the relative error in radius that encloses 
95% of the results. The radius values of the grid points have been upwardly adjusted to incorporate 
conservative expectations of signal suppression at each period. 



We also evaluate the overall detection probability for putative transiting planets. Given the 
cadence and coverage of the £^POX/ observations, we determine the number of in-transit points we 
expect for a given radius, period, and phase (where the phase is evaluated from to 1 periods, in 
increments of 30 minutes). Wc then evaluate the expected significance of the detection, assuming 
a boxcar-shaped transit at the depth of (i?p/i?^)^, and the noise of the actual GJ 436 time series. 
At each phase and period individually, we scale down Rp to incorporate the signal suppression at 
that ephemeris. We use the improvement in the over the null hypothesis to define a positive 
detection, after we have removed the best candidate transit signal (presented in Section 4.1). If 
we do not first remove this signal, then we are a priori defining a "detectable" signal to be any 
signal more prominent than the best candidate signal, and we would be unable to evaluate this 
signal's authenticity. We set our detection limit at an improvement in over the null hypothesis 
of 250. This level is set by the results of the Monte Carlo analysis; we observe a improvement 
of 250 compared to the null hypothesis only at radii and periods at which we correctly recover 
the injected period (with the exception again of the single signal which presented only one transit 
event). Although we use the ^ ^ detection statistic, we do not assign a percent significance 
to a Ax^ of 250 because of the presence of red noise in the time series. We then determine the 
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fraction of phases for which the signal would be significant enough for EPOXI to detect, which 
is shown in Figure [H For planets larger than about 2.0 times the size of Earth, we would detect 
the planet with 95% probability for all periods less than 8.5 days, and we would detect the planet 
with nearly 100% probability at periods less than 6 days — this includes the 2:1 resonance with the 
known Neptune-mass planet. We would also have a good chance (80% certainty) of detecting a 1.5 

planet with a period corresponding to the 2:1 resonance or shorter. Our sensitivity to 1 
planets is limited to periods of 1.5 days or less for greater than 50% confidence. Although planets 
with longer periods are much less likely to transit, we are still sensitive to planets larger than 2.0 
R@ with 80% certainty at a period of 15 days, and 70% certainty at a period of 20 days. 




Period [Days] 



Fig. 4. — Detection probability versus period for planets ranging in size from 1 to 2 i?0. The 
detection criteria is set by the percentage of phases at a given period for which the number of 
points observed in transit produces a improvement of 250, compared to the null hypothesis. We 
assume a boxcar-shaped transit at the depth of (Rp/Ri,)'^, where Rp at each period and phase 
includes the expected level of signal suppression. 



3.3. Dynamical Constraints on Second Planets 

Over the past few years, many researchers have pointed out that the excellent precision of 
transit data can reveal second planets. Such a detection would rely on the deviation of the known 
planet from a Keplerian orbit; in particular, the orbit and light curve would no longer be perfectly 
periodic. We may sort such effects into several categories. 

First, on orbital timescales the two (or more) planets exchange small amounts of angular mo- 
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mentum and energy as they approach and recede from n iutual conjunctio i is. This effect is imprinte d 
on fluctuations of the transit-to-transit orbital period (jAgol et al.ll2005l : iHolman &: Murravl l2005l ) . 
which might be seen in a short string of consecutive transits, like those we are reporting. Each of 
the 8 EPOCh transit times of GJ 436b is plotted in Figure O and given in Table El We consider 
variation around a constant period of > 100 s to be ruled out by our data. We introduce a small 
second planet in a series of numerical integrations spanning the transits that OX/ measured, to 
see what mass planet would have been detectable. For each period of a hypothetical planet of mass 
0.01 Mjup, we chose a circular orbit, coplanar with the known planet, and a grid of initial orbital 
phases. We recorded the time of minimum separations between the star and planet b, representing 
the EPOXI mid-transit times, and found their standard deviation from a linear ephemeris. We 
examined the minimum and the median (on the grid of phases) of those standard deviations. The 
minimum was generally considerably higher than the limits from other effects, quoted below. The 
median, however, was smaller than 100 s in orbits close to planet b. Even where it was not, as the 
short-term perturbations are linear in the mass of the perturber, we scale the value of the computed 
perturbations and the mass to 100 s and a "typical" mass that could have been seen in our data. 
We plot that typical mass in Figure [6] as a thin line. 
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Fig. 5. — Transit timing residuals from the new ephemeris. The solid diamonds are the previously 
published times from sources listed in the text; the hollow diamonds are the EPOCh times. The 
propagated error in the new period is overplotted as dashed lines. 



Second, on timescales of several tens of orbital periods (for the mass of GJ 436b), orbit-orbit 
resonances between planets can cause the periods to oscillate, the eccentricities to fluctuate, and 
the apses to precess. Our data alone does not have the requisite time span to sample all of such 
oscillations, although resonant effects do cause the short-term constraints of Figure [6] to appear 
jagged, as a function of period. In combination with light curves taken over several seasons, the data 
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gathered by EPOXI allow very deep constraints on planets in resonance. We defer this analysis 
to future work, as sampling of the relevant timescales is no t yet complete using available data. 
Similar analyses have alr eady been perf ormed on TrES-lb (Isteffen &: Ago! 2005), HP 189733b 
(jMiller-Rjcdetall boOSbl ) . HD 20 9458b JAgol &: SteffenlboOTi : iMiller-Ricci et al.1 boOSal ) . CoRoT- 
Exo-lb (jBeanlbood ). and TrES-3b (joibson et ahlbood ). The only considerable difference in the case 
of GJ 436b is that, because of its small mass, its response to resonant perturbations is proportionally 
bigger, so the current non-detections of perturbations allow surprisingly low-mass planets to be 
ruled out. For instance, once the 6 month timescale is well sampled by transit times with precisions 
of 7 seconds, then the data co uld detect or rule out a > 0.05MEarth planet librating with large 
amplitude in the 2:1 resonance (jPont et al.ll2009l ). 



Third, on timescales of thousands of orbital periods, the planets torque each other's orbits into 
different orientations and eccentricities, exchanging angular momentum, but not energy. These 
secular effects therefore do not change the semimajor axes and periods of the planets, but several 
recent papers have shown that these changes would still manifest themselves in transit data. For 
some effects, an eccentric planet is required and a large impact parameter is desirable, so GJ 436b 
(with an impact parameter of 0.85) provides a suitable laboratory to search for them. The remainder 
of this section evaluates these constraints in detail. 

Apsidal motion causes the longitude of transit to move either towards the periastron or towards 
the apastron of the orbit, depending on a; and the sign of Co (here uj is the angle between the 
ascending node on the plane of the sky and the periastron of the orbit). The precession period is 
expected to be considerably longer than the observational baseline, Pprec = 27r/a; ^> Tobs, so we 
expect to detect at most a linear change in oj. One consequence of precession is that the observed 
period of transits generally differs from the true (sidereal) orbital period. However, because the 
true period is known only from radial velocity observations, which are several orders of magnitude 
less precise than transit data for this purpose, in practice this method is not constraining. On 
the other hand, the period of transits would differ from the period of occultations, so multi-epoch 
secondary eclipse observations should be able to place a tight constraint on Co from this technique 
(jMiralda-Escudell2002l : lHevl &: Gladmanll2007l ). Currently, the technique with the highest precision 
is to use transit durations to constrain Co, as the duration changes linearly in time depending on 
the true planet-star separation at transit. This dependence comes in two forms: (1) the impact 
parameter is directly proportional to this separation and (2) the tangential velocit y is inversely 
proportional to this separation, according to Kepler's second law. It has been shown (|Kopallll959l ) 
that these two effects cancel, for linear changes in cj, at an impact parameter of 6 = l/\/2. At the 
impact parameter o f GJ 436b, the depen dence of the impact parameter is more important, as has 
been emphasized bv lShporer et al.l (|2009l ). Both effects together cause the duration to change at a 
rate: 

duj dT /dtl + esmoj 1 — 6^ 



e cos OJ 1 — 26^ 



(2) 



(Pal fc Kocsii 



2008 



Jordan &: Bakosll2008l ). where T is the duration between the times when the 
projected centers of the planet and star are separated by Ri,. We combined the transit duration of 
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our combined light curve, considered as a "single epoch" measurement of extremely high quality, 
with the transit durations reported over the last t wo years. Here and b elow, we take the values 
from the compilation and homogeneous analysis of ICoughhn et alJ (l2008l )FI We measure dT/dt = 
(1.9 ± 3.7) X 10"'^ min/day, so infer Co = (—8 it 15) x 10~^ deg/day by equation [51 Taking |a;| < 
4.5 X 10"'^ deg/day to be a conservative uppe r limit, we now consider what constraints this puts on 



additiona l planets in the GJ 43 6 system. As lHevl &: GladmanI (j2007l ) pointed out, the asymptotic 
formulae (jMiralda-Escudell2002l ) for apsidal rate due to second planets underestimates their strength 
by a factor of > 2 when the perturbing planet is within a factor of 2 in semi- major axis of the 
transiting planet. Thus we use the full formula for apsidal motion c aused by a perturbing planet 



on a circular orbit ( Hevl &: Gladman 



2007 



Murrav Dermottlll99g ): 



duj vr Mc 
'dt ~ IPb'Ml 



< ab 
> ab 



where 



^3/2 



x) 



27r 



COS tpdip 



(3) 



(4) 



iQ (1 + 2x cos i/j + x)^/2 

is the Laplace coefficient, which we compute numerically. We may invert this formula to obtain a 
constraint on the mass of a second planet, as a function of its semimajor axis, given the previously 
quoted constraint on w; this is plotted in Figure [6l Perturbing planets on currently eccentric orbits 
can torque planet b at somewhat different rates, depending on the relative apsidal orientation, an 
effect we neglect here. 

Another secular effect is nodal precession, which changes both the node and the inclination 
of the orbit. The former is measurable in principle by comparing Rossiter-McL aughlin (R-M) 
obser vations spanning many years, but GJ 436 is a slow rotator {Pmt ~ 48 days; iDemorv et al 



20071 ) and thus not a favorable R-M target. The latter, however, is exquisitely measurable due 
to the large impact p arameter of GJ 436b. In fact, detect ions of an inclination change have 
already been suggested (jRibas et al ] l2008l :l ICoughlin et al which would imply a several Earth- 

mass plai iet in a nearby orb i t, and the specific orbits and masses of all owable planets have been 



2008 



Demorv et al. 



200S 



Bean &: SeifahrtI 120081 ). We measured di/dt 



debated (jAlonso et al. 

(—6 lb 31) X 10~^ deg/day by producing a weighted l inear e ph emeris to our incl in ation and the 



ones 



presented in 



Couehhn et al.l (2008), 



Gillon et al 



(|2007al M. iBean et all (120081 ). IPeming et al 



(j2007l ). IShporer et al.l (j2009l ). and lAlonso et al.l (j2008l ). We find a a conservative upper limit of 
\di/dt\ < 9 X 10"^ deg/day. The shape of the light curve is much more sensitive to inclination 
change than to a psidal motion. Next, we con vert this to a mass constraint on second planets. The 
secular theory of iMurrav &: DermottI (|l999l . §7.2), gives the rate of inclination change due to a 



perturbing planet, as the nodes of the planets precess: 



Coughlin et al.l (|2008l ') analyzed the data in two ways; we use the results assuming fixed stellar and planetary 



radii for all the data sets. 



where duo/dt is given by Equation [3] and AOsky is the ascending node of the second planet relative 
to the ascending node of the known planet, measured clockwise on the plane of the sky. (This latter 
quantity is the sky-projected mutual inclination, similar to the sky-projected spin-orbit angle from 
the R-M literature.) Equation [S] assumes small eccentricities and mutual inclination. We have 
tested it against numerical integrations for the known planet and a second planet that would be 
marginally detectable and is far from resonance (Mc = 35M^, Pc = 2.818-P;,), and find it to be 
accurate to 20% over the timespan of the data if the perturbing planet begins on a circular orbit and 
is inclined less than 20° from the orbit of planet b. The resulting mass limit is plotted in Figure[6]for 
l^f^skyl ^ 10°. Because of the lightcurve's sensitivity to inclination change, only a small AQsky is 
needed for an inclination change over a certain angle to be as detectable as apsidal motion through 
the same angle. For AJlgky = 0.37°, the inclination constraint has the same magnitude as the 
apsidal motion constraint. According to Equation [5] it also has the same form. 

We note a degeneracy between the effects of apsidal motion and nodal precession. The light 
curve is most sensitive to both via transit duration change, so wit h a judiciously chosen orbit for 



a second planet, the two effects mostly cancel at the present time (|Miralda-Escudell2002l ). In such 
a case, the precession of the planet in its plane causes the impact parameter to increase at the 
same rate the precession of the planet out of its plane causes the impact parameter to decrease. 
This happens for AOsky = —0.37°, so a much larger Mc is possible for that particular orbit. The 
degeneracy could be broken by the shape of the light curve, as the changing velocity accompanying 
apsidal motion is not canceled, so the change in the ingress/egress duration and the depth of the 
transit due to the limb-darkened star may be observable. However, in practice the radial velocity 
constraint provides a stronger limit for such orbits. 

We have also plotted as a gray a rea in Figure [6] t he orbits near GJ 436b that are not guaranteed 



to be stable by Hill's criterion (e.g., lGladmanlll993l ). We took account of the eccentricity of b, so 
the stability criterion is approximately 



'^outer 'dinner 9 , r ^7/1 / "^^inncr ~l" -^outer \ ^/"^ /„\ 

> Woefo + 5.761 1 . (6 



There is only a narrow region that both satisfies the stability criterion and is ruled out by 
these secular constraints more strongly than by the radial velocity measurements. However, these 
dynamical constraints will certainly tighten in the future with both longer observing baselines and 
the measurement of the occultation period, to be compared with the transit period. The secular 
limits given above will improve approximately linearly in time, until 1 radian of the secular cycle 
has elapsed. This timescale may be computed using Equation [3] with max(M6,Mc) substituted 
for M;,, and is typically decades. In contrast, most relevant short-term and resonant timescales 
are shorter than the several-year time span of data collected so far, so apart from sampling all 
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the timescales, those constraints have saturated. Thus, they will only improve considerably with 
improved precision of transit data. 

In the case of GJ 436, excellent radial velocity data was in hand before transits were ever 
discovered, so it is not surprising that radial velocity constraints are stronger than the transit 
timing constraints for a wide range of periods of perturbers. However, soon exquisite transit data 
will be available from the Kepler mission, with little or no useful constraints from radial velocity. 
Therefore, the techniques illustrated here for GJ 436 are expected to be very useful in those cases. 

Finally, there are still several potentially confounding effects which are well below the sen- 
sitivity of the data. The apsidal line may precess not by a second planet, but by either the 
post- Newtonian relativistic correction, by a t idal bulge raised on the planet (jRagozzine Sz Woli 



20091 ). or by rotational oblatenes s of the star 

magnitude 6.5 x 10"^ d eg/day (|Pal Kocs id l2008l: I Jordan k Bakod I2008l 1. 1.0 x 10 



Miralda-Escudg 120021') . These effects have expected 

^ deg/day 



( Ragozzine fc Wolj[2009l l. and 3 x 10"^ deg/dav (jMiralda-Escudelbood )" respectively, all far below 
the current reach of the data. Also, nodal precession may occur and be observable if the stellar 
spin axis and the orbit normal do not lie in the same angle as projected on the plane of the sky. 
Such precession is driven by rotational oblateness, so its timescale is similar to apsidal motion by 
the same mechanism, approximately 3 x 10^^ deg/day. 

We note also that a detailed se arch for Trojan astero ids to the known planet G J 436b using the 
transit timing method proposed bv lFord &: Gaudil (|2006l ) could place interesting limits on Trojans 
in this system, although we do not specifically address this question in this work. However, Trojan 
bodies could also produce a photometric transit, so our constraints on additional transiting bodies 
could be applied to the detection of Trojans as well. These Trojan asteroids would have a period 
equal to GJ 436b (trailing or preceding the hot Neptune planet at the fixed phase set by the L4 
and L5 Lagrange points), so without conducting a detailed analysis specific to Trojan asteroids, we 
conservatively would have been able to detect Trojans that produced a transit depth equivalent to 
a 1.5 planet with nearly 100% certainty from the detection probabilities shown in Figure HI 



3.4. Rotation Period of GJ 436 

As previously stated, we fit a sinusoidal modulation in the GJ 436 light curve with a polynomial 
and remove it in order to conduct the search for additional transits. Here we consider the astrophys- 
ical significance of this signal, which we attribute to changes in the apparent brightness of the star 
due to star spots. The sinusoidal m odulation has a pea k amplitude near 0.5 millimagnitudes and is 



shown in the top panel of Figure[71 I Butler et al.l ( 20041 ) give an upper limit on V sin i of 3 km s~^; 
if we assume rigid body rotation and a stellar radius from this work of = 0.437 it 0.016 Rq, 
we infer a lower limit on the rotation period of ~ 7.4 days. We investigate the periodicity of the 
EPOCh GJ 436 observations by binning the time series in two hour increments and creating a 
Lomb-Scargle periodogram of the binned observations. We find only one peak with false alarm 
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Fig. 6. — Limits in the mass - period plane for a hypothetical companion planet GJ 436c. The 
secular limits are given due to a lack of observed transit inclination and duration change, assuming 
the companion is currently on a circular orbit. A coplanar second planet must have a mass below 
the long-dashed line, due to our upper limit on apsidal motion of b. A non-coplanar second planet 
must have M2 below the thick solid line for |Ar2sky| ^ 10°, with the limit scaled up or down in pro- 
portion with the actual value of |Ar2sky|. Such planets could be stable; according to Hill's criterion 
(|Gladmanlll993l ) second planets on initially circular orbits are certainly stable outside of the gray 
area. The approximate radial velocity limit is given as a critical sem i-amplitude K (short-d ashed 
line) that would have been detectable according to the simulations of lBean &: SeifahrtI (|2008l ). For 
perturbers producing short-term fluctuations that would have been detectable in EPOXI alone, a 
typical value of the detectable mass is given as a thin solid line; in contrast to the other lines, this 
should not be viewed as an upper limit (see text). 



probability significantly less than 0.01% at Prot=9.01 days; the time series is shown phased to this 



period in the bottom panel of [7 
limit set by V sin i. However, 



The best-fit rotatio nal period of 9.01 days is longer than the lower 



Demorv et al.l (|2007l ) observed flux variations of GJ 436 nearly 20 



times larger with amplitude of 1% which were best fit by a stellar rotation period « 48 days. It 
is therefore possible that the baseline of the EPOXI observations is not long enough to resolve the 
full rotation period of the star, and the modulation we observe is due to multiple starspots passing 
into view over a period shorter than the rotation period of the star. 
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Fig. 7. — Top panel: GJ 436 out-of-transit and out-of-eclipse EPOCh observations, binned in two 
hour increments. Middle panel: Lomb-Scargle periodogram of binned observations, with false 
alarm probability of 0.1% overplotted. The most significant peak corresponds to a period of 9.01 
days. Bottom panel: Time series phased to most significant period of 9.01 days. 
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4. Discussion 



4.1. Best Candidate Transit Signal 

We find no transit signals for GJ 436c at the significance limit set by Ax^=250. We present 
our best candidate here, which corresponds to a planet with radius 1.04 i?® and period 8.42 days. 
The improvement in the over the null hypothesis is 170. We investigate the possibility of signal 
suppression by masking these suspected transit points from the points used to create the 2D spline, 
and recreating the spline surface. We find that the candidate transit depth is unaltered, indicating 
that these points are well-sampled on the CCD. In Figure [HI we show the phased and binned best 
candidate signal, as well as the two observed candidate transit events separately. The gap in the 
middle of the first event is due to the star wandering off the bottom edge of the CCD for a period 
of about an hour. 
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Fig. 8. — Top panel: Best transit candidate, shown phased (black points) with best-fit model 
overplotted in red. The green curve directly above is the phased light curve binned by a factor 
of 30. Bottom panels: The two observed candidate transit events which contribute to the best 
candidate signal, with best-fit model overplotted in red. The transit event at left contributes an 
improvement over the null of 49, and the right transit contributes an improvement in of 121. 
Since we place our significance limit at Ax^ = 250, we cannot claim a positive detection. 



We also perform an identical transit search on a flipped version of the time series (we subtract 
1 from the normalized time series and multiply by -1), since we expect red noise fluctuations to 
introduce both positive and negative imposter transit signals. We find that the best solution to 
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the inverted time series produces an improvement over the null hypothesis of Ax^ 



106. 



4.2. Radius constraints 



From the results of our Monte Carlo analysis and phase coverage analysis, we rule out transiting 
planets orbiting interior to GJ 436b larger than 1.5 Earth radii with 95% confidence. We would 
expect to detect such a planet with 95% certainty, and recover its radius to within 15% with 95% 
confidence (referencing Figures H] and [3l respectively). For planets exterior to the orbit of GJ 436, 
we can no longer assume that additional planets will transit the host star — since GJ 436b itself has 
an inclination of i = 86.60 it 0.17°, and the host star has a radius of R^, = 0.437 it 0.016 (from 
this work), we would anticipate coplanar planets to transit only with periods less than about 3.4 
days. However, the orbital inclinatio n of Earth differs from the orbital inclinations of the gas giants 
by up to 2.5° in the case of Saturn (ICoxll2000l ). If the putative GJ 436c had an inclination which 
was 2° closer to edge-on than GJ 436b, we should be able to detect transits out to nearly 13 days. 
For all periods less than about 8.5 days, we would have detected the transit of a planet 2.0 i?® or 
larger with 95% probability if the planet produced a transit (shown in Figured]). This includes the 
2:1 resonance with GJ 436b at about 5.3 days, although NICMOS observations of the transit times 
of GJ 436b, with timing variations less than a few seconds, disfavor planets in the 2:1 resonance- 
a planet in the 2:1 resona nce with GJ 436 b with a mass as small as 0.01 should produce at 
least 7 second variations (jPont et al.ll2009l ). For periods less than about 20 days, we would have 
detected the transit of a planet of 2.0 or larger with 70% probability. 



4.3. Mass constraints 



We use the theoret ical mass-radius relationships for super Earths and hot Neptunes calculated 
bv lFortnev et al.l (|2007l ) to place approximate mass constraints, given our radius constraints derived 
from the search for additional transiting planets. We rule out planets larger than 1.5 -R® interior 
to G J 436b — using the analytic formulae given in iFortnev et al.l (j2007l ) , such a planet would have 
a mass around 0.8 if it were pure ice, 2.9 if it were pure rock, and 4.6 if it were pure 
iron. At semi-major axes larger than that of GJ 436b, since additional planets may not transit, 
we are unable to set firm upper limits from the lack of transits. If such a planet with a period 
less than 8.5 days did transit, we would be sensitive with 95% confidence to radii as small as 2.0 
R^- even at a period of 20 days, although the transit probability is much less likely, we would still 
detect a planet this size with 70% certainty. We can therefore rule out transiting planets at periods 
less than 8.5 days with masses greater than 2.3 assuming a pure ice composition, 9.6 
assuming a pure roc k composition. Alth ough a pure iron planet with a radius of 2.0 R^ would have 
a mass of 63.5 (jFortnev et al.l 120071 ). this con iposition is perhap s unrealistic, even assuming a 
formation history with mantle-stripping co llisions (IMarcus et al ] l20ld ). The maximum mass of a 2.0 
Rq planet, using the relations derived by Marcus et al. ( 20ld 'l for maximum collisional stripping. 
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would be closer to 20 M®. Planets with periods less than 7 days have been ruled out in the G J 436 



syste m by radial velocity constraints down to about 8 with 3a confidence (jBean &: Seifahrt 



20081 ). Our limits on the presence of pure rock planets are therefore complementary with previous, 



stronger constraints. 

We find that the EPOCh dynamical constraints on additional planets with periods from 0.5 to 
9 days rule out coplanar secular perturbers as small as 10 and non-coplanar secular perturbers 
as small as 1 M^, as shown in Figure [6l These dynamical constraints are not as strong as current 
radial velocity constraints, except in orbits very close to that of GJ 436b. However, we anticipate 
that dynamical analyses similar to those presented in this work will prove useful to the community 
in cases of planets with masses below current radial velocity detect ability, such as those that Kepler 
will find. 



4.4. Eccentricity of GJ 436b 



The eccentricity of GJ 436b has been attributed to two possible mechanisms. First, the residual 
eccentricity can be attributed to e xcitation from th e dyna mical interactions of GJ 436b with an as- 
yet undetected additional planet. iBean &: SeifahrtI (j2008l ) tested this second hypothesis by finding 
how well the radial velocity data could be improved by adding perturbers to the system and evolving 
the system forward by numerically integrating the Newtonian equations of motion. Their analysis 
ruled out perturbers greater than 8 at periods less than about 11 days (semi-major axes less 
than 0.075 AU) with high confidence. They presented radial velocity solutions that improved the 
fit by up to 4(7 at smaller masses with periods between 4 and 11 days. We rule out rocky transiting 
bodies down to 9.6 with periods less than 8.5 days with 95% confidence in the GJ 436 system. 
However, th is doesn't prec l ude t he possibility of an additional planet at these periods which does 
not transit. I B at v gin et al.l (|2009l ) compiled a list of possible dynamically stable secular perturbers 
which are consistent with the transit times, radial velocities, and observed eccentricity of GJ 436b. 
These proposed additional planets all have periods greater than 16 days, however. EPOXI would 
be sensitive to a transiting planet larger than 2.0 i?® with close to 80% probability at 15 days (see 
Figure SD . 

The second possible explanation for the eccentricity of GJ 436b is a tidal Q parameter that is 
much larger than that of the ice giants in our solar system (thereby increasing the circularization 
time to greater than the age of the syste m) — such a Q would need to be 1-2 orders of magnitude 
larger than that measured for Neptune (jBatvgin et al.ll2009l : iBanfield &: Murravlll992l ). However, 
the tidal Q for Jupiter may be as high as 2 x 10^ with the assumption that Jupiter and lo are 
orbiting in a steady-state configurati on; and may be eve n higher i f that assumpt i on is false (and 
tides on lo are currently dominant) (| Jackson et al.ll2008l ). In fact, lOgilvie &: LinI (|2004l ) find that 
hot Jupiters may typically have Q values near 5 x 10^, so a value of 10^'^ for GJ 436b for Q/k2, 
proposed bv I Jackson et al.l (l2008|'l. mav no t be unreas onable (the Love num ber k2 is typically near 
0.5 for Solar System gas giants: lBursalll992l ). However. iBatvgin et al.l (j2009l ) suggest that the actual 
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Q for GJ 436b must be higher still. 



A definitive explanation for the eccentri city of GJ 436b is so f ar undetermined, but a resolution 
to this question is observationally tractable. iBatvgin et al.l 1)20091 ) provide a thorough discussion of 
follow-up observations that could measure the signal of a secular perturber to GJ 436b: radial ve- 
locity measurements are sensitive enough at their current level to resolve the periodogram signature 
of such a per turber, and a long b aseline of transit times at the level of At < 10 s could also confirm 
its presence (jBatvgin et al.ll2009l ). If these methods find no signal corresponding to a perturbing 
body, then the tidal Q of GJ 436b may ind eed be much h i gher t han those for the ice giants in our 
solar system. The explanation proposed by I Jackson et al.l (|2008l ) may alternatively be correct- the 
Q values for solar system giant planets are current underestimations of the true Q value for these 
planets, in which case the tidal Q necessary to explain GJ 436b may not be inconsistent with that 
of Neptune. 
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Table 3. EPOCh Derived Stellar and Planetary Parameters for GJ 436 



Parameter 



Value 



Light curve Parameters 
P (days) 

Transit Times (BJD) 



R*/a 
i (deg) 



2.64389579 ± 0.00000080 
2, 454, 592.76134 ± 0.00022 
2, 454, 595.40458 ± 0.00026 
2, 454, 598.04822 ± 0.00021 
2, 454, 600.69271 ± 0.00045 
2, 454, 603.33657 ± 0.00017 
2, 454, 605.98079 ± 0.00035 
2, 454, 608.62326 ± 0.00026 
2,454, 611.27015 ±0.00040 
0.08142 ± 0.00085 
0.0707 ± 0.0025 
86.60 ±0.17 



Stellar and Planet Parameters'' 

(Rq) 0.437 ±0.016 

Rp (R®) 3.880 ±0.147 



^Using the stellar mass from I Torre 



